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ABSTRACT 


The  finite  element  method  is  applied  to  The 
two-dimensional,  inviscid,  compressible  radial 
equilibrium  equation  for  axial  compressors. 
Isoparametric  elements  are  used  along  with  three-point 
Gaussian  integration  for  stiffness  matrix  evaluation. 
The  radial  equilibrium  equation  is  put  into 
quasi-harmonic  form  for  stream  function  formulation 
and  results  are  presented  using  an  isentropic  flow 
assumption.  Axial  velocity  profiles  at  rotor  and 
stator  blade  edges  are  compared  with  published 
performance  data  of  the  NASA  Task-1  stage  transonic 
compressor  and  with  numerical  finite  element  results 
of  Hirsch  and  Warzee. 
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I.   INTRODUCTION 


A.   PROBLEM  STATEMENT  AND  OBJECTIVE 


The  prediction  of  meridional  flows  within  turbo  machines, 
be  they  compressors  or  turbines,  is  a  difficult  but 
important  part  of  the  design  process.  The  difficulty  arises 
from  the  presence  of  three-  dimensional  and  viscous  effects 
within  all  turboraachines  and  the  importance  arises  from  the 
necessity  to  design  accurately  and  efficiently. 

To  simplify  the  problem  of  viscous,  three-dimensional 
analysis,  Wu  [ Ref .  1  ]  showed  that  this  complicated  flow  may 
be  analyzed  by  solving  two  interrelated  flows:  one  being  the 
blade-to-blade  flow  describing  the  flow  between  rotating 
blades  and  the  other  being  the  meridional  through  flow  which 
describes  the  radial  equilibrium.  These  flows  are  depicted 
in  Fig  1.  In  addition,  an  inviscid  and  axi  symmetric 
assumption  is  made  in  the  through-flow  thereby  simplifying 
the  flow  to  a  two-dimensional,  axi  symmetric,  inviscid,  and 
compressible  analysis. 


Three  methods  may  be  found  in  current  reports  regarding 
the  solution  of  the  radial  equilibrium  equation.  The  first 
two  are  the  streamline  curvature  method  [Ref. 2, 3, and  4]  and 
the  matrix  method  [Ref. 5  and  6]  which  is  oasically  a  finite 
difference  technique.  The  third,  a  relatively  new  method,  is 
the  finite  element  method.  As  shown  by  Hirsch  and  Warzee 
[Ref. 7]',  the  solution  of  the  radial  equilibrium  equation  by 
the   finite   element   method   is   achieved   by  arranging  the 


equation  for  the  stream  function  in  quasi-harmonic  form. 

Due  to  the  excellent  results  reported  in  Raf.7  and  to 
further  the  research  effort  for  finite  element  techniques  in 
fluid  flow  problems,  the  purpose  of  this  thesis  is  two  fold. 
Firstly,  the  goal  was  to  formulate  a  computer  program  for 
solution  of  the  radial  equilibrium  equation  paralleling  the 
steps  as  presented  by  Hirsch  and  Warzee.  Secondly,  after 
suitable  verification  of  computer  results  with  those  of 
Hirsch  and  Warzse,  the  goal  was  to  compare  computer 
predicted  flows  with  measured  performance  data  of  the  Naval 
Post  Graduate  School's  transonic  compressor. 


The  purpose  of  this  paper  is  to  present  a  report  on  the 
results  obtained  thus  far.  In  Section  II,  the  derivation  of 
the  radial  equilibrium  equation  is  presented  followed  by  the 
application  of  the  finite  element  method  to  this  equation. 
Section  III  describes  the  computer  program  in  some  detail. 
Section  IV  contains  selected  test  cases  which  wer=  used  in 
program  testing  and  checking.  In  Section  V,  conclusions  are 
presented  along  with  recommendations  for  further  study  and 
work  on  the  project.  The  appendices  contain  the  program 
listing  along  with  a  sample  test  case  for  reference  by  the 
user.  In  addition,  a  list  of  references  is  contained  for 
further  reading  on  the  subject  of  this  paper. 
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Figure    1    -      MERIDIONAL  AND   BLADE- TO- BLADE    PLANES 


II.   THEORY 


A.   THE  DERIVATION  OF  THE  RADIAL  EQUILIBRIUM  EQUATION 


The  following  discussion  is  taken  from  Ref.  7  with 
slight  changes  in  notation.  The  basic  turbomachine  geometry 
to  be  analyzed  is  depicted  in  Fig  2.  Although  the  machine 
noted  is  one  stage  of  a  compressor,  a  similar  analysis  to 
the  one  that  follows  may  be  applied  to  other  machines  such 
as  axial  turbines  and  mixed-flow  machines. 

One  begins  with  the  Euler  equation  assuming  the  viscous 
forces  to  be  negligible. 

— -> 


^  +  (v-V)  /-  Vp/f        (ii. a.  1) 


The  continuity  equation,  assuming  unsteady  flow  is, 


^    +  V(fV,)=  O  (II.  A.  2) 


The   First   Law   of   thermodynamics   in   a   fluid   field 
becomes, 


Tvs-  vh-  vp/f  (II.&.3, 
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Figure  2  -   TURBOMACHINE  GEOMETRY 
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Substituting  equa 
leads  to  the  Crocco  equation, 


tioQ  (II. A. 3)   into   equation   (II. A. D 


(II.  A.  4) 


where  B  is  the  total  enthalpy. 

steady   and   adiabatic   flow,   the   energy 
Assuming   a   steaay   auu 

equation  becomes  simply, 

(v-v^o  <II-A-5) 

which   shows  that  along  a  streamline  in  a  stationary  system, 
the  total  enthalpy  is  constant. 

IB  a  relative  system  such  a_s  the  case  In  a  coto. :  blade 
ro«,  the  total  relative  velocity,-,  can  he  expressed  rn  -he 
following  form, 

where  5   is  the  constant  angular   velocity   and  S      is  the 

a    «-p  -t-ho  relative  system, 
constant  peripheral  speed  of  the  reiati 

^i-i-^p  svstem  becomes, 
Now,  the  Crocco  equation  m  a  relative  system 

-7  7\         -r-      nfi     W1  _  l^l2  ]   (II. A. 7) 

,TT  .  5^  for  the  stationary  system,  the 
Parallel  to  equation  (II.i.5)  tor  x  ,rPiativ*) 

energy  equation,  assuming  steady   and   adiabatic   (relative) 
flow  in  a  relative  system,  become 
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s 


(w-  7)Hr  = 


(II.  A.  8) 


e  =o 


where  H  is  the  relative  total  enthalpy  expresed  as  follows, 


i  oz 


H„3  ^    w*   _  ^»g 

i     ~2 


(II.  A.  9) 


From  the  following  velocity  diagram, 


Vm  >  wm 
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equation  (II. A. 9)  may  be  arranged  as  follows. 


Since, 


w^+  wze  =  wz  «Ym  HV*  -U) 


(II. A. 10) 


then  , 


w%  Y^  ^v,v  -ZUVo  +U: 


(II.  A.  1  1) 


and, 


w*-  v%  ^-  zu\/* 


(II. A.  12) 


Substituting    equation    (II. A. 12)     into      equation       (II. A. 9) 
leads    to    the    following   relation, 


«*-  h  +  t  -UY*  =    H-W, 


e- 


(II. A.  13) 


Equation   (II. A. 8)   shows   that   "    is   constant   along    a 
streamline  in  a  relative  system. 


Upon  integrating  equation  (II. A. 8)  between  the  rotor 
inlet  and  outlet,  the  Suler  equation  for  turbonachines  is 
found , 


AH 


OU~i 


\N 


A  (U-Y*J 


0U.I 


,N 


(II. A.  14) 
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It  may  be  shown  [Ref.9]  that  by  circumf erentially 
averaging  equation  (II. A. 1),  and  under  the  axi  symmetric 
flow  assumption  the  following  relation  is  valid, 


-Y<(V^)-  T-  V5  -  VH  -t^   +  Fj 


(II. A. 15) 


where  Fi  is  the  body  force  of  the  blades  acting  on  the  fluid 
and  all  variables  are  mean  values  along  the  direction  of  the 
circumference.  Hence  equation  (II. A. 15)  is  an  approximation 
for  axi  symmetric  flow.  As  a  final  note  on  equation 
(II. A.  15),  since  the  viscous  forces  were  neglected  in 
equation  (II. A. 1) ,  there  must  be  a  force  introducing  the 
entropy  variations  along  the  blade.  This  force  is 
proportional  to  the  pressure  loss  coefficient  and  is  labeled 
Fj  ,  the  dissipative  force.  F^  produces  work  which  in  turn 
produces  entropy  production  radially  along  the  blade.  Under 
the  axi  symmetric  assumption,  entropy  varies  axiaily  and 
radially  only  and  is  assumed  to  be  proportional  to  the 
pressure  loss  coefficients  [ Ref .  2  and  8]. 

•  Due  to  boundary  conditions  imposed  on  the  problem  and 
the  axi  symmetric  assumption,  cylindrical 
coordina tes,  {r,&, zj  , will  be  used  in  all  subsequent  analysis. 
Therefore,  equation  (II. A. 15),  in  cylindrical  coordinates 
and  with  axial  symmetry  is  as  follows, 


^K)-V^v<-^yJ^-T^  -^-f-     ."-Li «» 


%  X   A,  ,  1    V£  X 


tt(to)+  7fi<?V.J»  F* 


<e^ 


(II.  A.  17) 
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Ve  t&k  -  k\)  -  %H)  -  H  -  T  ^  -  F?     (i i .  a .  1  a, 


It  is  important  to  note  here  that  under  the  axi 
symmetric  assumptions,  equation  (II. A. 15)  reduces  to  the 
following, 


V'V- 


o 


(II. A. 19) 


Likewise  in  a  relative  system  (rotor)  ,  the  axi  symmetric 
assumption  leads  to  the  following, 


— ">  ->» 


t  =  o 


(II. A. 20) 


Equation  (II. A. 16)  describes  the  meridional  through  flow 
radial  equilibrium  equation  for  the  finite  element  method. 
Since  one  is  concerned  with  the  meridional  plane,  the 
following  derivative  expression  is  taken  from  Fig  3. 


Yv^^wv"  Vjt  ££    "*■  Vj  J^_ 


(II. A. 21) 


Therefore    equation     (II. A. 17)     reduces    to, 


x 


^e-   ^%,  (Mo) 


(II. A. 22) 


which  reveals  that  in  a  duct  where  there  are  no   blades   and 
therefore   no   blade   forces,   angular   momentum  is  constant 


1b 


along  a  streamline.  In  that  case, 


^(KVpJ-0 


(II. A. 23) 


As    shown   in   Ref.9,   the   circumf erentially   averaged 
continuity  equation  is  the  following, 


Ji(<\™*)  +  k{fW-° 


(II. A. 24) 


where  b  is  the  blockage  factor  defined  by  Hirsch  and  Warzee 
as  the  tangential  area  reduction  due  to  the  thickness  of  the 
blade. 


k-i-4- 


(II. A. 25) 


where  t  is  blade  thickness  and  s  is  blade  spacing. 
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R 


Figure  3  -   MERIDIONAL  PLANE 


IS 


One  further  step  in  the  formulation  of  the  radial 
equilibrium  equation  for  solution  by  the  finite  element 
method  involves  introducinq  the  stream  function.  In 
cylindrical  coordinates,  the  stream  functions  are  defined  as 
follows, 


V  -  -J-   ^ 


(II. A. 26) 


(II. A. 27) 


Substitutinq   these  expressions  into  equation  (II. A. 16), 
the  equation  becomes, 


(II. A. 28) 


The  riqht  hand  side  of  equation  (II. A. 28)  is  applicable 
to  the  absolute  flows  in  the  stator  and  duct  reqions.  For 
relative  flows  such  as  those  in  the  rotor,  the  riqht  hand 
side  is  modified  by  replacinq  the  total  enthalpy, H,  by  the 
relative  total  enthalpy, H,  and  the  quantity  Ve/R  is 
replaced  by  He/R . 

As  a  last  assumption  in  the  formulation  of  the  qoverninq 
relation  for  the  meridional  throuqh  flow  radial  squilibrium 
equation,  both  the  radial  component  of  the  body  force,  Fj, 
and  the  radial  component  of  the  dissipative  force, F^  ,  are 
neqiected.  This  assumption,  [Ref.1,8]  does  not  hamper  the 
accuracy  of  the  results  for  conditions  at  desiqn  speed. 
Even   thouqh   published  compressor  performance  data  used  for 
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the  test  case  in  this  thesis  was  obtained  at  0.5  design 
speed,  these  force  terms  were  also  neglected  in  the  computer 
program.  As  will  be  shown  later,  this  assumption  could 
possibly  have  had  adverse  effects  on  the  predicted  axial 
velocity  profiles  at  the  rotor  hub  and  tip  regions. 


The  final  representation  of  the  meridional  radial 
equilibrium  equation  to  be  solved  by  the  finite  element 
method  is  as  follows, 


where , 


A.  f,   JL+  "\         i.    (,    JL+ 


3it(k3*  )+  M^]  +  j  =0       l"-4-29' 


k-    ^£b 


(II. A. 30) 


and 


(II.  A.  31) 


B.   THE  FINITE  ELEMENT  METHOD  APPLIED  TO  THE  RADIAL 
EQUILIBRIUM  EQUATION 


In  order  to  formulate  equations  (II. A. 29)  through 
(II. A.  31)  in  matrix  form  for  solution  by  the  finite  element 
method,  one  must  apply  a  weighted  residual  technique  to  the 
equations  for  numerical  solution.  The  weighted  residual 
method  used  here  is  the  Salericin's  Method.  The  following 
discussion  is  taken  from  Ref.  7  with  only  slight  changes  in 
notation . 

Rewriting  equation  (II. A. 29)  and  dividing  through  by  R, 
one  has, 
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UtUkl^ 


(II.  B.  1) 


=  0 


where  this  equation  represents  the  flow  in  the  volume, V. 

The  boundary  condition   for   this   partial   differential 
equation,  after  dividing  through  by  R,  is, 


-fc  k^  +*,  OM.)    =o  (II-B-2) 


where  this  equation  solves  the  flow  on  the   close!   boundary 
of  the  volume, or, S. 


3y  applying  the  weighted  residual  process  to  equations 
(II.B.1)  and  (II. B. 2)  and  using  an  arbitrary  weighting 
function,  w(r,z),  one  has 


J  Vj  (<ii)  rVol  Av  +   j  w.aj  r5oB  d  £'    =  o 


(II.  B.3) 


where   r    and   r„     are   the   volume  and  surface  residuals 
respectively,  or, 


^--\\^)^M)  + 


-  o 


(II. B. 4) 


y 


5ur  -   £ 


Krt   W.^^o)] 


(II. B. 5) 


If      the      solution      to    equation    (II.B.1)     was    exact,     both    r 
and    r„    _    would    be    equal    to    zero. 


Vol 


In   order   to   clarify   the  boundary  condition,  equation 
(II. B. 2) ,  one  may  analyze  the  equation  as  follows. 
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On  the  surface 


e, S, where  y   is  specified, 


<M, 


(II. B. 6) 


and, 


<   - 


i  ~^  oo 


(II.  B.  7) 


Similarly,  on  the  surface,  where  )yi  " u » S  ,  where 


4 1    =0 


(II. B. 3) 


£.n&=o,     £,  u£  =  S* 


(II. B. 9) 


Due  to  the  axi  symmetric  assumption,  the  final  equation 
will  not  involve  dV  and  ds  but  the  intersection  of  dV  and  dS 
with  the  meridional  plane.  Therefore,  one  must  transform  the 
volume  integral, dV,  to  a  surface  integral  and  the  surface 
integral, dS,  to  a  line  integral. 


and, 


Hence,  let, 

qJI  =  intersection  of  dV  and  meridional  plane 
qC  =  intersection  of  dS  and  meridional  plane 

<j£r  2Tr£dC 

With   this   transformation,   one   may   rewrite   equation 
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(II.  B.  3)     as    follows, 


-wot^UV^SVila.cliit 


*+ 


where    on   the   contour,    C,  V=T©« 


-  o 


(II. B. 10) 


One  must   now   integrate   the   first   term   in   equation 
(II.B.10)  by  parts  to  obtain  the  following, 


«£_ 


(II. B.  1  1) 


Inspecting  the    first    term    in    equation     (II.B.11) /    one   may 
use    the    following    integral    theorem    to    simplify    further 


Ji^c/Jl    r      fcjmdC 


(II. 3. 12) 


Jl 


Rewriting      equation    (II.B.11)     ,gives, 


/■  J  J  (II. 3.1 


3) 


4-UlK^ZndC  =  0 


Finally,    since 


^  ^ 


tf 


^Y\ 


ae 


2% 


(II. B.  14) 


equation  (II.B.13)  reduces  to  the  following, 
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-a  '44 


(II.  B.  15) 


One  now  has  the  final  equation  in  the  form  for  use  by 
the  weighted  residual  method  using  any  arbitrary  weighting 
function,  W  (r,z)  .  As  noted  previously,  the  Galerkin's 
Method  will  be  used  here  which  implies  that  the  weighting 
functions  are  the  same  functions  used  in  approximating  the 
stream  function,  t  • 


Before  applying  the  finite  element  method,  one  must 
discretize  the  continuum  and  then  approximate  the  unknown 
function  ,Y  ,  by  a  set  of  polynomials.  For  this  particular 
problem,  eight-noded  iso-parametric  elements  were  chosen  for 
discretization,  see  Fig  4,  and  the  following  approximating 
functions  were  used. 


i 


(II.B.16) 


where, 

M"  l^7!  J  =  s^aPe  functions 

%     =  value  off  at  the  node 

f  =  value  of  4^  at  an7  arbitrary 
location  within  the  element.  The  shape  functions , N- ,  used 
here  are  defined  by  the  following  relations  as  shown  in 
Ref . 10, 
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Mi(^)  =    iO-f^'  +  Vl\0  (II. 3. 17) 


where  the  following    coordinate   transformations   are    used, 

S 

(II  .B.  18) 

U.    (&    Vi\a 
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**    I.^Hi^ 
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Figure   4    -       ISOPARAMETRIC   QUADRILATERAL  ELEMENT 


2b 


At  this  point,  one  is  ready  to  apply  the  Galerkin  Method 
to  equation  (II.B.15)  by  substituting  equation  (II.B.16)  for 
the  unknown  function^  /  an^  N-  for  the  weight  function,  W, 
which  yields, 


K^i^V  ^I^Va  -  \i*  ** =  0(ii-b-i9: 


-*l 


c-i 


t~( 


-a. 


This         integration        yields        the      following      system      of 
equations    which    is    solved    for    the    unknown    nodalH^r 


where, 


[ 


r, 


l^u   X\z  "'  K\v\ 

kv\i  •  •  •   k*m 


n  H 


(II.B.20) 


(II. 3. 21) 


and, 


I-   [fM* 


(II.B.22) 


In  addition  since  both  the  'stiffness  matrix', K,  and  the 
right  hand  side  vector,  [F],  are  functions  of  V  ,  the  system 
as  defined  by  equations  (II. 3. 20)  through  (II. 3. 21)  must  be 
solved  iteratively. 

At  this  point  one  has  the  total  finite  element 
formulation  of  the  radial  equilibrium  equation  as  defined  by 
equations  (II.B.19)  and  (II. 3. 20).  The  problems  which 
remain  to  be  clarified  are  basically  two  fold.   Firstly  one 
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must  evaluate  the  integrals  in  equations  (II.B.19)  and 
(II. B. 2)  by  numerical  methods,  and  secondly,  the  solution 
procedure  for  the  non-linearity  must  be  formulated.  In  Part 
C,  both  of  these  final  steps  are  presented. 


C.   NUMERICAL  INTEGRATION  OF  STIFFNESS  MATRIX  AND  SOLUTION 
PROCEDURE 


1.   Numerical  integration  of  the  stiffness  matrix 

As  noted  in  Section  II.  B,  evaluation  of  equation 
(II.B.21)  must  be  performed  numerically.  In  addition,  one 
realizes  that  the  derivative  expressions  enclosed  within  the 
interval  must  be  evaluated  by  a  coordinate  transformation. 
This  is  done  in  the  following  way, 

Since, 


i_ 


(II. C.  1) 


then , 


^ 
^ 


^M;  &    ^-Nc  Xr 


&   * 


t 


}r  *-< 


(II.  C.  2) 
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and  in  matrix  form, 
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(II. C. 3) 


Furthermore,  defining  the  Jacobian  matrix  as, 


r 


J  , 


A* 

±r 

1 

i 

it 

>• 

*l 

*l 

(II. C. 4) 


then   by   dividing  both  sides  of  equation  (II. C. 3)  by  J,  one 
has  the  following  transf ocmation, 


-[ 


vi 


J 


f^i 


2 


(II. C. 5) 


V 


In    addition,    it    has    been    shown    [Ref.9]    that 


J*clr    =   Wl<^cl>| 


(II. C. 6) 


Now,   with   equations   (II. C. 5)    and    (II. C. 6),    equation 
(II.B.21)  becomes  the  following, 


29 


Mi 


(II. C. 7) 


Equation  (II. C. 7)  is  best  integrated  using  the 
Gauss-Legendre  integration  method  since  it  is  of  the 
following  form, 


»  i 


*ij.  j  [e^i^i 


or  finally,  [Ref.10], 


^rjMH.V 


(II.  C.  8) 


(II. C. 9) 


where  A. and  B. are  coefficients  (Fig   5)   for   both   two  and 
three  point  Gaussian  Quadrature. 

At  this  point,  one  has  the  tools  to  calculate  all 
the  elements  of  the  stiffness  matrix.  In  like  manner,  the 
right  hand  side  vector, f,  is  calcualted  by  numerical 
integration . 


30 


NUMBER  OF 

GAUSSIAN 

POINTS 

+  A- 

+  B- 

2 

0.57735  02691 

1.00000  00000 

3 

0.77459  66692 
0.00000  00000 

0.55555  55555 
0.88888  88888 

Figure  5  -   GAUSS  IAN  INTEGRATION  POINTS 
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2  •   Solution  procedure 

The  following  is  a  synopsis  of  the  basic  solution 
process.  Specific  details  concerning  equations  and  methods 
of  computer  coding  are  covered  in  the  proceeding  section. 
The  proceeding  is  meant  to  give  the  reader  a  preview  of  the 
solution  process. 

a.  Discretization 

Initially  the  machine  under  analysis  is 
discretized  into  eight-node  iso-parame t ric  elements.  The 
axial  calculation  stations  are  placed  arbitrarily  in  the 
duct  regions  and  along  blade  edges  and  centers  for  the  rotor 
and  stator  as  shown  in  Fig  6.  At  this  point  the  system 
topology  and  nodal  coordinates  are  specified. 

b.  Initialization 

To  begin  the  iteration  process,  one  must  assume 
an  initial  internal  stream  function,  velocity,  and  density 
distribution.  In  the  program,  the  initial  internal  stream 
function  was  assumed  to  be  that  of  the  outer  boundary 
throughout  while  the  velocity  and  density  distribution  was 
assumed  to  be  that  of  the  inlet. 
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Figure  6  -   COMPRESSOR  DISCRETIZATION 
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c.   Calculation  of  thermodynamic  variables 

Before  calculating  the  right-hand  side  vector, f 
,  one  must  obtain  distributions  of  angular  momentum, 
enthalpy,  and  entropy.  This  is  done  by  first  calculating 
the  thermodynamic  variables  at  the  inlet  axial  station  from 
the  given  inlet  conditions.  In  order  to  proceed  axially 
through  the  machine  to  calculate  the  nodal  angular  momentum, 
enthalpy,  and  entropy,  the  following  three  equations  derived 
in  Section  III  are  used. 


H  -  CpT  =  constant  along  a  stator  streamline 


fyr  Cftr 


(ujx) 


=  constant  along  a  rotor  streamline 


X  V©  =  constant  along  a  duct  streamline 


An  example  of  this  calculation  procedure  for  the  iuct  region 
is  shown  graphically  in  Fig  7. 
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Figure    7    -      DUCT  ELEMENT 
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In   this  figure,  the  angular  momentum  at  point  X 
is  equal  to  the  angular  momentum  at  point  Y.  Wore  formally, 


^e)r-l3^,\)H\-  {r\l,\ 


(II.  C.  2.  1) 


Since   the   previous   axial    station's   thermodynamic 
variables   are    known,    one    must    now    find    the    values    of     ^       and 
v^     at      point      Y.    This    is    done    iteratively    in    the    following 
way.    Since 


(II. C.  2.2) 


and  along  the  left  side  of  the  element, 


\ 


=  -1 


(II.  C.  2.  3) 


then  equation  (II.C.2.2)  may  be  solved  for  Yj  by  a  suitable 
iteration  method.  As  will  be  shown  in  the  next  section,  a 
half-interval  method  was  used  to  obtain  the  unknown  Y[ 
Once  VI  is  known,  then  equation  (II. C.  2.1)  is  solved  for 
the  angular  momentum  at  point  X.  The  rotor  and  stator  are 
handled  in  a  similar  fashion.  In  addition,  the  rotor  and 
stator  deviate  the  flow  creating  a  three-dimensional  flow 
field  between  the  blades  in  the  respective  blade  row.  Low 
speed  cascade  correlation  data  [Ref.13]  was  used  to 
calculate  the  effective  turning  angles  in  the  rotor  and 
stator.  These  effects  are  calculated  beforehand  with  known 
mass  flow  rate  and  unifora  axial  velocity  assumptions  at  the 
rotor  inlet.  The  results  of  these  calculations  are  part  of 
the  input  data  routine  in  the  form  of  relative  and  absolute 
flow  angles  at  the  rotor  nodes  and  absolute  flow   angles   at 
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the   stator   nodes.   This   will  be  shown  more  exactly  in  the 
next  section. 


d.  Calculate  matrices 

At  this  point,  the  right  hand  side  vector, f,  and 
the  stiffness  matrix  ,K,  are  calculated. 

e.  Solve  system  of  equations 

The  system  of  equations  as  shown  in  equation 
(II.B.20)  is  solved  for  the  nodal  stream  function. 

f.  Perform  relaxation  iteration 

Due  to  the  strong  non-linear  properties  of  the 
system  of  of  equations,  the  following  iterative  scheme  is 
necessar  y. 


(II.C2.tO 


where  c(  is  the  under  relaxation  factor.   As  will  be  shown  in 

Section  III,   this   scheme   is   performed   only  in   certain 

regions   of   the   machine   and  in  addition  after  a  specified 
number  of  iterations. 

g.   Update  velocity  and  density  profiles 

Using  the  current  nodal  distribution  of  the 
stream  function,  axial  and  radial  nodal  velocity  components 
are  calculated  along  with  a  new  nodal  density   distribution. 
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Again,   this  calculation  procedure  will  be  shown  in  the  next 
section. 


h.   Test  for  convergence  of  Y 


Stream  function  convergence  criteria  is  now 
tested  and  will  determine  if  further  iterations  are 
necessary.  The  solution  is  said  to  converge  if  the  following 
equation  holds  for  all  nodes. 


%   -  X 


<v 


r>  +  i 


<  £ 


(II. C. 2. 5) 


where  €   is  a  designated  requirement  for  convergence. 

i.   Summary 

In  summary,  the   eight   steps   involved   in   the 
solution  are  noted  below; 

(1)  Discretize  the  continuum. 

< 2)  Assume  an  initial  stream  function,  velocity, 
and  density  solution. 

(3)  Calculate   the  nodal  thermodynamic  variables 
from  the  given  inlet  conditions. 

(4) Form   the  right  hand  side  vector,  f(r,z),  and 
the  stiffness  matrix,  K. 


(5)  Solve  the  system  of  equations,  given  by,  [K] 
=  [F]  for  a  new  stream  function  distribution. 


(6)  Perform  relaxation  iteration  if  required. 

(7)  Calculate  new  nodal  velocity  and  density 
distributions  from  the  current  stream  function  solution. 

(8)  Test  the  solution  for  convergencs,  and  if 
required,  repeat  steps  (3)  through  (8)  using  the  current 
nodal  stream  function  values. 

This  concludes  the  solution  description  and  now 
one  is  ready  to  more  completely  understand  the  computer 
program  which  assembles  the  preceding  eight  steps. 


39 


III.   THE  PROGRAM 


A.   OVERALL  FLOWCHART  AND  DESCRIPTION 


The  overall  flowchart  of  the  program  is  depicted  in  Fig 
8.  Those  blocks  denoted  by  the  letter  '  S'  are  subroutines, 
while  the  remaining  calculations  are  an  integral  part  of  the 
main  program. 

After  proper  dimensioning  of  all  arrays  and  subseguent 
initialization,  the  input  data  are  read  and  then  printed. 
This  not  only  presents  a  physical  picture  of  the  problem  but 
also  serves  as  a  cross  check  to  the  user  for  correct  data 
insertion.  In  addition,  a  subroutine  is  availabla  to  obtain 
a  computer  drawn  plot  of  the  mesh  (Fig  6)  and  is  a  further 
check  on  proper  data  input. 

At  this  point  ail  the  necessary  variables  have  been 
stored  and  the  iteration  counter  for  stream  function 
convergence  is  set.  tfith  the  current  nodal  values  of  t  and 
the  given  inlet  thermodynamic  conditions,  the  thermodynamic 
variables  throughout  the  machine  are  calculated.  From  the 
calculated  values  of  enthalpy  and  angular  momentum, 
(isentropic  flow  is  assuaed) ,  the  right-hand  side  vector  is 
calculated  followed  by  the  stiffness  matrix  calculation 
(eguation  II.B.21). 

The  system  of  equations  (equation  (II.B.20))  is  now 
solved  for  the  new  nodal  stream  function  distribution.  It 
is  here  where  for   all   iterations   but   the   first   that   a 
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relaxation  factor  is  applied  as  noted  previously  in  equation 
(II. C. 2-4).  The  reasoning  behind  not  applying  the 
relaxation  scheme  to  the  value  of  nodal  ^  after  the  first 
iteration  is  the  fact  that  the  first  iteration  produced  a 
close  approximation  to  the  correct  stream  function 
distribution.  With  this  close  approximation  to  the  stream 
function  came  a  velocity  and  density  distribution  which  in 
turn  was  near  the  correct  solution.  It  was  found  that  if 
the  first  iteration  was  relaxed,  the  second  iteration  became 
unstable  since  in  fact  the  velocities  and  densities  were 
themselves  farther  from  the  true  values  than  wsre  assumed 
initially. 

After  testing  the  nodal  stream  function  for  convergence 
by  use  of  equation  (II.C.2.5),  the  calculation  process  is 
either  repeated  or  ceased  by  virtue  of  convergence  or 
limiting  the  number  of  iterations. 

As  stated  previously,  low  speed  cascade  correlation  data 
[Ref.13]  were  used  to  calculate  turning  angles  in  the  blade 
regions.  These  angles  were  assumed  constant  throughout  the 
solution  and  not  refined  after  subsequent  iterations. 
Further  work  on  the  computer  program  could  entail  an 
additional  computational  routine  which  would  calculate  the 
new  turning  angles  after  each  iteration.  A  sample 
calculation  of  rotor  turning  angles  is  shown  in  Appendix   D. 

In  the  following  sections  the  program  structure  is 
examined  in  more  detail. 
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DIMENSION  ARRAYS 


INITIALIZE  ARRAYS 


READ   IN   DATA 


SET  ITERATION  COUNTER 


CALCULATE  THERMODYNAMIC 
VARIABLES  FROM  GIVEN  f 


CALCULATE  RIGHT  HAND 
SIDE  VECTOR  FROM  GIVEN 
THERMODYNAMIC  VARIABLES 


CALCULATE  STIFFNESS  MATRIX 


SOLVE  SYSTEM  OF  EQUATIONS  FOR  f 


PERFORM  RELAXATION 
ITERATION  IF  NOT  1st  ITERATION 


CALCULATE  NEW  VELOCITY  AND  DENSITY 
DISTRIBUTION  FROM  NEW  STREAM  FUNCTION 


PRINT  FINITE 
ELEMENT  RESULTS 


Figure   3  -   PROGRAM  FLOWCHART   (STOP) 
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B.   TH2  MAIN  PROGRAM 


1 •   The  input  routine 

The  following  is  a  description  of  the  input  data 
required  by  the  program.  The  data  are  arranged  into  twelve 
categories  described  in  the  following  manner. 

a.  category  1 

Problem  identification. 

b.  category  2 

Number  of  nodes  and  number  of  elements. 

c.  category  3 

Node    numbers,    nodal   coordinates   and   nodal 
blockage  factor. 

d.  category  4 
System  topology. 

e.  category  5 
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Element  type;  duct,  rotor,  or  stator. 


f.   category  6 


Absolute  flow  angles  for  rotor  and  stator  nodes, 


g.   category  7 


Relative  flow  angles  for  rotor  nodes, 


h.   category  8 


Inlet  thermodynamic  quantities. 


i.   category  9 


Physical  constants  for  fluid  under   observation, 


j.   category  10 


First  estimate  of  internal  stream  function 


k.   category  11 


Node   numbers    and    specified    nodal    stream 


function . 


1.   category  12 
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Node  numbers  where  the  right  hand  side,  f  (r,z), 
is  to  be  calculated. 

Before  describing  in  detail  the  format  to  be 
followed  for  data  insertion,  it  is  important  to  note  the 
following  assumptions. 


(1)  Uniform  flow   conditions   at   inlet   and 


outlet. 


(2)  Uniform  flow  conditions  at  rotor  inlet 
for  calculation  of  appropriate  turning  angles.  This 
assumption  is  necessary  to  calculate  the  values  of  rotor  and 
stator  flow  angles. 


continue. 


detail. 


With   this   in   mind   ,   the    discussion    will 


The  following  describes  each   category   in   more 


Category  1 : 


Format:  (20A4) 


Number  of  cards:  1 


Procedure:  Enter  the  title  of  the  problem  in 


columns  1-20. 


Category  2: 

Format:  (2110) 

Number   of   carls:   Equal   to   the  number  of 
nodes  in  the  system. 
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Procedure:  Enter  the  number  of  nodes  in 
columns  1-10,  and  the  number  of  elements  in  11-20.  Both 
integers  must  be  right  justified. 

Category  3: 

Format:  (I10,3F10.0) 

Number  of  cards:  Equal  to  the  number  of 
nodes  in  the  system. 

Procedure:  Each  card  contains  the  node 
number  followed  by  the  Z  coordinate,  R  coordinate,  and  nodal 
blockage  factor.  The  coordinates  are  in  dimensions  of 
inches. 

Category  4: 

Format:  (915) 


Number  of  cards: 
elements  in  the  system. 


Eaual   to   the   number   of 


Procedure:  Each  card  contains  nine  integers 
right  justified  in  columns  5,  10,  15,  etc.,  through  45.  The 
first  integer  is  the  element  number  followed  by  the  eight 
nodes  associated  with  that  element.  It  is  important  to  note 
that  the  nodes  are  read  in  starting  with  the  upper  right 
hand  node  and  proceeding  in  a  counterclockwise  fashion 
around  the  element. 

Category  5: 

Format:  (2110) 

Number   of   cards:   Equal   to   the  number  of 
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elements. 

Procedure:    Enter   the   element   number   in 
columns  1-10,   followed   by   the   integer   '1'   (duct) ,   '2' 
(rotor),  or  '3'  (stator)  describing  the  element  as  either  in 
a  duct,  rotor,  or  stator  region. 

Category  6: 

Format:  (6 X , A4 ,1  10 , F10  .  0) 

Number  of  cards:  Equal  to  the  number  of 
rotor  and  stator  nodes  plus  one  'STOP'  card. 

Procedure:  Enter  the  node  number  (right 
justified)  in  columns  11-20  followed  by  the  value  of  the 
associated  absolute  flow  angle  in  radians  in  columns  21-30. 
The  last  card  in  this  category  is  a  'STOP'  card  entered  in 
columns  7-10. 

Category  7: 

Format:  (6 X , A 4  , 1  10,F 1 0  .  0) 

Number  of  cards:  Equal  to  the  number  of 
rotor  nodes  plus  one  'SID?'  card. 

Procedure:   Enter   the   node   number  (right 

justified)   in   columns   11-20   followed  by  the  value  of  the 

associated  relative  flow  angle  in  radians  in  columns  21-30. 
The  last  card  in  this  category  is  a  'STOP'  card. 

Category  8: 

Format:   (7F10.0)  ,  (F10.0) 
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Number  of  cards:  2 

Procedure:  Enter  the  following  quantities  in 
the  prescribed  order  and  with  the  noted  dimensions. 

First  card 

Mass  flow  rate:  (lbm/sec) 

Inlet  axial  velocity:  (ft/sec) 

Outlet  axial  velocity:  (ft/sec) 

Inlet  total  density:  (lbm/ft  ) 


Inlet  static  density:  (lbm/ft  ) 


Inlet  total  pressure:  (lbf/in  ) 


Inlet  total  temperature:  (*R) 


Secoad  card 


Speed:  (RPS) 


Category  9: 


Format:  (3F10.0) 


Number  of  cards:  1 


Procedure:  Enter  the  following  quantities  in 
the  prescribed  order. 

Gas  constant:  (f t-lbf/lbm- *R) 
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Ratio  of  specific  heats 


Constant 


pressure     specific      heat 


(BTU/lbm-  R) 


Category  10: 

Format:  (F10.0) 

Number  of  cards:  1 

Procedure:  Enter  the  first  estimate  of  the 
internal  stream  functiion  to  be  used  in  the  first  iteration. 

Category  11: 

Format:  (5  X , A4 , 1 10,F 1 0  .  0) 

Number  of  cards:  Equal  to  the  number  of 
nodes  having  a  specified  value  of  the  stream  function  plus  a 
•STOP'  card. 

Procedure:  This  set  of  cards  allows  the 
stream  function  boundary  conditions  to  be  read  in.  A 
typical  card  contains  an  integer,  right  justified  in  columns 
11-20  ,  which  is  the  node  number,  followed  by  the  value  of 
the  specified  stream  function  in  columns  21-30.  The  last 
card  is  a  'STOP'  card. 

Category  12: 

Format:  (6X,A4,I10) 

Number  of  cards:  Equal  to  the  number  of 
nodes  where  the  right  hand  side  is  to  be  specified. 
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Procedure:   Enter   the   node   number,  right 

justified   in  columns  11-20,  where  the  right  hand  side  is  to 

be  calculated.   Again,  the  last  card  in  this  category  is   a 
•STOP1  card. 

After  all  the  data  has  been  read  by  the  program, 
the  input  data  is  printed  and  the  mesh  is  plotted  for 
verification  by  the  user.  The  sample  format  is  shown  in 
Appendix  C. 

This  concludes   the   input   routine.    The   next 
section   describes   the  calculation  of  the  stiffness  matrix, 

K. 

2.       Stiffness,   matrix    evaluation 


As  sho 


wn  previously  in  Section  II.C.1,  the  following 


equation  describes  each 


term  in  the  eight  by  eight  elemental 


matrix. 


^-\\^fi\W\\[^[^\ 


(III.B.2. 1) 


in   addition,    •*•    is   defined    in    the    following    way   in    order   to 
numerically    integrate   the    equation. 

i 

k  =   — 5 S  '  _  (III.B.2. 2) 

L-KI      »  «-  =  « 

j.~-i    k  i  Ar  if  a  n  p  factor  taken  as 
wherp  b  is  defined  as  the  elemental  blockage  tacto 

u+    „„*<><*    of   tr«   particular   element 
an  average  over  the  eight  nodes  or       v 

c       a      -an«-ouadraturs   integration 
and   <^)   are   the   defined   ,auss  Quaira 

points . 
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Figure  9  -   STIFFNESS  MATRIX  EVALUATION 
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Fig  9  depicts  the  flowchart  for  both  elemental 
stiffness  matrix  evaluation  and  the  assemblage  into  the 
system  stiffness  matrix.  More  specifically,  the  figure 
shows  a  three-point  Gaussian  Quadrature  scheme  but  can  be 
changed  to  a  two-point  scheme  by  simply  integrating  four 
times  instead  of  nine  as  shown. 

The  actual  coding  of  the  stiffness  matrix  evaluation 
and  assemblage  may  be  found  in  lines  STR03510  through 
STR0U770  in  the  computer  program. 

3 .   Solution  of  systems  of  e^uat ions 

At  this  point,  the  system  of  equations  are  modified 
for  the  boundary  conditions  and  solved  for  the  nodal  stream 
function  values.  An  equation  solving  routine,  DSI.1Q, 
available  in  the  system  library  was  used  for  this  purpose. 
It  was  found  that  no  comparable  savings  was  realised  by 
using  a  banded  equation  solver. 

4«   Iteration  schemes 

As  noted  previously  in  Section  II. C,  a  relaxation 
scheme  is  necessary  for  convergence  to  a  solution. 

Two  distinct  differences  with  regard  to  the 
iteration  method  were  noted  from  that  of  Ref.7.  Firstly,  it 
was  found  that  relaxation  was  necessary  only  in  the  rotor 
and  stator  elements  and  also  in  the  duct  region  between  the 
rotor  outlet  and  stator  inlet.  Secondly,  due  to  the  extreme 
non  linearity  in  the  rotor-stator  areas,  a  switch  was 
required  which  changed  the  sign  of  <(  in  equation  (II. C.  2.4) 
as  required  for  stability  of  convergence.   Clarification   of 
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this  change  follows:  It  was  found  that  during  the  initial 
three  or  four  iterations,  the  stream  function  values  of  the 
rotor-stator  nodes  sometimes  exceeded  the  value  of  the  upper 
boundary.  Due  to  an  abscence  of  sources  within  the  domain 
of  solution,  this  occurrence  was  incompatible  with  the 
boundary  conditions.  At  this  point,  it  was  necessary  to 
make  cL  negative  in  equation  (II.C.2.4).  During  subsequent 
iterations,  as  the  solution  converged,  the  rDtor-stator 
regions  became  stable  and  the  sign  of  c(.  was  returned  to  its 
positive  value.  This  iteration  proved  to  stabilize  the 
solution  with  respect  to  stream  function  values  and 
velocities. 

The  iteration  procedure  is  coded  in  the  computer 
program  from  lines  STR05070  through  STR05200. 

5-   The  output  routine 

Once  convergence  is  obtained  or  the  number  of 
iterations  have  reached  the  limit  imposed  by  the  user,  the 
results  are  displayed.  A.  sample  output  is  shown  in  the 
Appendix.  In  addition,  the  units  of  all  dependent  variables 
are  the  same  as  those  noted  in  the  input  routine. 


C.   THE  SUBROUTINES 


The  following  describes  each  of  the  six  subroutines  in 
the  computer  program.  Each  subsection  contains  a  list  of 
calling  arguements  and  for  subroutines  ECAL,  SLINE,  and  VEL, 
a  basic  flowchart.  In  addition,  for  those  subroutines  whose 
mathematical  theory  was  not  presented  in  Section  III,  a 
brief  treatment  is  also  given. 
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1 •   Subroutine  shape 

This  subroutine  calculates  the  shape  functions 
(equation  (II.B.17))  at  the  values  of  ^  and  Y)  as  requested 
in  the  arguement  list  below. 

SUBROUTINE  SHAPE  (E,Z,SF) 
E  =  value  of  $  (input) 
Z  =  value  of  yi  (input) 
SF  =  eight  by  one  vector  of  the  eight  shape  functions. 

2 •   Subroutine  Jacob 

JACOB  calculates  the  Jacobian  matrix  as  defined  in 
equat  ion  (II.  C.  1.  4)  for  the  value  of  %  ,Y|  denoted  in  the 
arguement  list. 

SUBROUTINE  JACOB  (E 1 , Z 1 , D, 2, RC$,ZC$ , RJAC) 

E1  =  value  of  yi  (input) 

Z1  =  value  of  §  (input) 

D  =  eight  by  one  vector  of  y«  (calculated) 

E  =  eight  by  one  vector  of  -j^r  (calculated) 

RC$  =  eight  by  one  vector  of  the  * r'  coordinates  of   the 
nodes  associated  with  the  element  (input) 

ZC$  =  eight  by  one  vector  of  the  'z'  coordinates  of   the 
nodes  associated  with  the  element  (input) 

RJAC  =  two  by  two  Jacobian  matrix  (output) 
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In  addition,  the  subroutine  assumes  that  the  vectors 
RC$  and  ZC$  contain  element  coordinates  arranged  in  a 
counter  clockwise  fashion  beginning  with  the  upper  right 
corner  node. 

3 .   Subroutine  si in e 

This  subroutine  calculates  the  thermodynamic 
variables  throughout  the  machine  given  the  inlet  conditions 
as  described  in  Section  II. C. 2.  The  calling  arguements  are 
defined  below. 

SUBROUTINE  SLINE(UINLET,RC,PSI,WRL,H,OVSL,VVEL,TVEL, 
NODE,NNODEI,CP,IT, KK,ALP,WG,TWEL,BS,HS) 

UINLET  =  Inlet  axial  velocity 

RC  =  Nodal  ' r'  coordinates  vector 

PSI  =  Nodal  stream  function  vector 

WRL  =  Nodal  angular  momentum  vector 

H  =  Nodal  total  enthalpy  vector 

UVEL  =  Nodal  axial  velocity  vector 

VVEL  =  Nodal  radial  velocity  vector 

TVEL  =  Nodal  absolute  tangential  velocity  vector 

NODE  =   Matrix   containing   nodes   associated   with   the 
element 

INLET  =  Vector  containing  node  nu-mbers  at  inlet  station 

NNODEI  =  Number  of  nodes  at  inlet  station 

CP  =  Specific  heat 

TT  =  Total  temperature  at  inlet 
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KK  =  Iteration  counter 

NTE  =  Element  type  vector 

ALP  =  Nodal  absolute  flow  angle  vector 

TWEL  =  Nodal  relative  tangential  velocity  vector 

BE  =  Nodal  relative  flow  angle  vector 

HS  =  Nodal  static  enthalpy  vector 

As  shown  in  Fig  10,  the  basic  calculation  procedure 
begins  with  calculating  the  reguired  energy  and  momentum 
values  at  the  inlet  station.  At  this  point,  beginning  with 
element  one,  the  element  type  is  interrogated  to  distinguish 
between  duct,  rotor,  and  stator  elements.  If  the  element  is 
in  a  duct  region,  then  the  streamline  intersections  for 
local  nodes  2,6,7,3  and  1  (Fig  7)  are  determined  along  with 
the  associated  values  of  energy  and  angular  momentum.  For 
the  rotor  and  stator  elements,  one  must  initially  find  the 
energy  and  momentum  values  at  local  nodes  3,4,5  (Fig  7)  due 
to  the  discontinuities  imposed  by  the  blade  edges.  Once 
these  calculations  are  performed,  then  the  process  for  the 
remaining  nodes  in  the  element  proceeds  in  a  similar  fashion 
to  the  duct  elements. 

After  all  the  elements  have  been  cycled  through,  the 
new  distributions  of  nodal  angular  momentum  and  energy  are 
returned  to  the  main  program  for  further  computations. 
Specifically,  these  values  will  be  used  by  the  next 
subroutine  ,  FCAL,  for  calculation  of  the  right  hand  side 
vector. 
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Figure    10    -      SUBROUTINE   SLINE 
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4.   Subroutine  fcal 


FCAL  calculates  the  right  hand  side  vector  as 
defined  by  equations  (II. A. 31)  and  (II.B.21).  Using  the 
identical  coordinate  transformations  for  numerical 
integration  as  described  in  Section  II. C,  the  final  equation 
to  be  coded  is  the  following, 


Mi 


-i-' 


(III.C.  4.  1) 


+  T(2,2)-MW-^^TJ1d\ 


where  isentropic  flow  is  assumed,  and, 
W-  =  angular  momentum 


H-  =  total  enthalpy 

v 


(III.C.  4.  2) 


The  arguement  list  is  defined  below.  In  addition, 
only  those  variables  in  the  list  which  have  not  been  defined 
previously  are  described. 

SUBROUTINE  FCAL (F , H , H , Z A , EA , U VEL, RC,ZC, AS L ,TV EL, NFS 
,NODE,NN,NE,NNFSP,TWEL,NrE) 

F  =  Right  hand  side  vector,  f  (r,z) 

W  =  vector  of  gaussian  quadrature  coefficients 

ZA  =  Vector  of  V  gaussian  quadrature  points 

EA  =  Vector  of  YL;  gaussian  quadrature  points 

NFS   =  Vector  containing  nodes  where  the  right  hand  side 
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is  to  be  specified 

NN  =  number  of  nodes 

NE  =  number  of  elements 

NNFSP   =   Number   of   nodes  where  the  right  hand  side  is 
specified 

Fig  11  depicts  the  basic  flowchart  for  the 
subroutine.  To  initialize  the  procedure,  one  begins  with 
the  first  node  (upper  right  hand  corner)  of  the  first 
element.  A  switch  is  then  applied  which  determines  if  the 
right  hand  side  is  to  be  calculated  at  the  node  or  if  a 
stream  function  value  has  been  specified.  This  information 
is  transferred  from  the  main  program  through  the  arguement 
list.  Once  the  node  is  allowed  through  the  switch,  then  the 
integration  process  is  started  at  the  first  integration 
point.  As  in  Section  III. A. 2,  the  flowchart  depicts  a 
three-point  Gauss  Quadrature  scheme.  After  the  integration 
has  been  completed,  a  switch  determines  if  all  the  local 
nodes  in  the  element  have  been  cycled  through  and  if  so, 
then  the  assembly  of  the  elemental  vector,  F$,  is  performed 
to  build  the  system  right  hand  side  vector  ,  F.  Finally, 
the  subroutine^determines  if  all  the  elements  have  been 
examined  in  order  to  signal  completion  of  the  right  hand 
side  vector.  At  this  point,  the  vector,  F(r,z),  is  returned 
to  the  main  program  for  problem  solution. 
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Figure  11  -   SUBROUTINE  FCAL 
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5 .   Subroutine  v  e  1 

This  subroutine  calculates  axial  and  radial 
velocities  and  also  densities  at  each  of  the  nodes  from  a 
known  stream  function  distribution.  As  noted  previously  in 
Section  II. C. 2,  both  velocity  and  density  profiles  are 
updated  after  obtaining  the  latest  value  of  nodal  stream 
function . 


The  velocity  calculation  proceeds   from   the   stream 
function  equations, 


V, .  -Li* 

v  -  -J- 4* 


(III.C.5.  1) 
(III.C.  5.  2) 


where  'b'  is  the  tangential  blockage  factor.   Since   r,  ,and 
are  of  the  following  form, 

t 


i 


^  i,  N;* 


(III.C.  5.  3) 


then  the  equation  for  the  axial  velocity,  V?  ,  becomes, 

1  r, 


V*« 


i-Mv 


S  ^H> 


kI^ M.I  KU  Li-.  >r  T- 


(III.C.  5.  4) 


Again,   since  the  shape  function,  N-  (t/i  )  ,  is  not  an 
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implicit  function  of  ' r '  and  'z',  one  must  use  equation 
(II.C.1.5)  to  obtain  the  proper  derivatives  for  computation 
of  equation  (III.C. 5. 4).  For  example,  from  equation 
(III.C.5)  , 


% 


Lr< 


:-i[r(l/)^+3-'(ui^] 


"i 


(III.C.  5.  5) 


At    this    point,    with     equation     (III.C. 5.5) 
substituted   into  equation  (III.C.5.4) ,  one  has  the  complete 
expression  for  the  axial  velocity  as  functions  of*  ,**].    One 
proceeds   similarly   for  expressinq  the  radial  velocity,  7  , 
in  terms  of  5  andYj  . 


In  order  to  calculate  the  nodal  density,  one  uses 
the  following  density  relation  for  flows  in  the  stator  and 
duct  regions. 


1 
•ft 


|_yu- 


(III.C. 5. 6) 


where 


£,  is  the  stagnation  density 


Since, 


v'*lv,,*v;Un-W'0 


(III.C.  5.7) 


then , 


H'-^^wA 


t1-! 


(III.C.  5.8) 
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Since  the  density  appears  or.  both  sides  of  the 
equation,  the  new  nodal  density  is  obtained  iteratively  at 
the  node. 

For  the  relative  flows  in  the  rotor,  the  following 
relation  for  static  density  is  used  [Ref.14], 


f 


VL1"^-? 


tt^k    _    {f-\\  W  *- £oN?l1  r-l  (III.C.5.9) 


<e 


Again,  the  solution  of  the  nodal  density  is  obtained 
in  an  iterative  fashion. 

In  the  following  arguement  list,  only  those 
variables  not  defined  in  the  previous  subroutine 
descriptions  are  noted. 

SUBROUTINE  VEL ( NE , NN, RC , NODE , G, EG , TT, RHOT , RHCN, ZC, 
PSI,RHO,B,UINLET,UVEL,VVEL,RHOSTA,NTE, ALP) 

G  =  Ratio  of  specific  heats 

RG  =  Gas  constant 

RHOT  =  Total  density  at  the  inlet 

RHON  =  Work  vector  which  contains  the  new  nodal   density 
distribution 

RHO  =  Nodal  static  density  vector 

B  =  nodal  blockage  factor  vector 

RHOSTA  =  Static  density  at  the  inlet  station 

The  basic  flowchart  for  SUBROUTINE  VEL  is  shown  in 
Fig  12.  Beginning  with  the  first  node  of  the  first  element, 
the   Jacobian   matrix   (equation  (II.C.1.4))  and  its  inverse 
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are  found.  At  this  point  the  partial  derivatives  with 
respect  to  •r*  and  ,z*  Df  the  shape  functions  are  found  as 
noted  in  equation  (III. C. 5. 5).  A  switch  then  allows  those 
nodes  not  at  the  inlet  station  to  pass  and  calculates  the 
new  density  and  velocities  at  the  nodes.  For  those  nodes  at 
the  inlet,  the  velocities  and  static  densities  are  retained 
at  the  given  inlet  conditions.  This  is  done  to  maintain 
boundary  condition  integrity  for  the  solution.  After 
cycling  through  all  elements,  the  subroutine  returns  the  new 
nodal  velocity  and  density  distributions  to  the  main 
program. 


o4 


YES 


V  =UINLET 
z 

V  =   0 
r 


ff.. 


in 


|  BEGIN  WITH   1st   ELEMENT,  11=7 


BEGIN  WITH   1st   NODE   OF   ELEMENT 
1=1 


11=11+1 


FIND  JACOB  IAN 


I  1=1+1 


FIND   INVERSE   OF   JACOB IAN 


FIND 


^r  4z 


FIND    If  Vp 


FIND      f 


n+1 


FIND  Vr,Vz 


YES 


RETURN 


NEXT   NODE 


[NEXT  ELEMENT 


Figure    12    -      SUBROUTINE    VEL 
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6  .   Subroutine  m_plo  t 

This  subroutine  utilizes  the  Calcomp  plotter  to 
depict  the  mesh  topology  of  the  machine  under  observation. 

SUBROUTINE  MPLOT (RC , ZC , NODE, NN , NE) 

This  completes  the  description  of  the  main  program 
and  associated  subroutines.  In  the  next  section,  a  test 
case  is  carried  through  from  data  input  to  final  results. 
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IV.   TEST  CASES  AND  RESULTS 


The  program  was  tested  by  using  published  performance 
data  [Ref.12]  of  the  NASA  Task-1  stage  transonic  compressor. 
The  compressor  was  discretized  into  twenty-eight  elements 
and  107  nodes  with  15  axial  calculation  stations  (Fig  6). 
The  speed  was  0.5  design  speed  with  a  mass  flow  of  107.6 
lbm/sec.  In  addition,  uniform  flow  was  assumed  both  at  the 
inlet  and  outlet  stations.  Turning  angles  for  the  rotor  and 
stator  were  pre  calculated  assuming  uniform  conditions  at 
the  rotor  inlet  and  using  NASA  SP-36  blade  correlation  data 
[Ref.13].  These  absolute  and  relative  flow  angles  were 
assumed  constant  throughout  the  iterative  procedure  as  they 
were  an  integral  part  of  the  input  data.  Tha  Appendix 
contains  a  listing  of  the  input  data  and  output  results  for 
the  NASA  Task-1  transonic  compressor  with  test  conditions 
noted.  To  compare  the  accuracy  of  the  predicted  flow  with 
actual  laboratory  observations,  computed  axial  velocity 
profiles  at  the  rotor  inlet,  rotor  outlet,  stator  inlet,  and 
stator  outlet  were  compared  with  experimental  results.  In 
addition  numerical  results  from  Ref.7  were  also  compared. 

Fig  13-16  show  the  computer  predictions  plotted  with  the 
experimental  values  and  the  numerical  solutions  obtained  by 
Hirsch  and  tfarzee.  The  profiles  shown  were  obtained  after 
ten  iterations  and  using  a  relaxation  factor  of  0.2.  The 
figures  show  that  the  best  overall  agreement  with 
experimental  data  occurred  in  the  stator  inlet  and  outlet. 
In  this  region  the  worst  error  was  17%  which  occurred  at  the 
stator  tip  inlet.  The  average  error  throughout  the  stator 
region  with  respect  to  experimental  data  was  6.6%. 
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The  rotor  hub  and  tip  outlet  area  exhibited 
instabilities  in  density  convergence  using  equation 
III. C. 5. 9.  Specifically,  the  density  solution  converged  to 
within  8%  at  the  rotor  outlet  tip  and  hub.  It  was  found 
that  by  not  allowing  the  nodal  density  at  these  nodes  to  go 
below  a  critical  value  of  0.06  lbm/cu  ft,  the  solution  for 
the  stream  function  converged.  By  allowing  the  nodal 
densities  at  the  rotor  outlet  tip  and  hub  to  go  below  this 
critical  value,  the  computed  velocities  at  these  nodes 
became  increasingly  large  and  the  arguement  within  the 
brackets  of  equation  III. C. 5. 9  became  less  than  one.  This 
prevented  continuation  of  the  iterations  for  the  stream 
function  solution.  In  addition,  the  rotor  tip  outlet 
exhibited  more  instability  than  the  rotor  hub  outlet.  The 
static  density  at  the  rotor  hub  outlet  oscillated  about  a 
value  of  0.062  lbm/cu  ft  while  the  rotor  tip  outlet  was 
constantly  driven  to  the  critical  value  of  0.06  lbm/cu  ft. 
One  method  attempted  to  alleviate  this  problem  was  the 
following.  Since  a  haLf-interval  iteration  routine  was 
used,  one  trial  run  involved  reversing  the  direction  of 
consecutive  guesses  when  the  density  iteration  did  not 
converge.  It  was  found  however,  that  after  three  to  four 
iterations  of  the  systea  of  equations,  the  static  densities 
at  the  rotor  outlet  tip  and  hub  were  again  driven  to  smaller 
and  smaller  values  which  led  to  instability  once  more.  The 
nodal  densities  converged  at  all  interior  points  of  the 
rotor  edge  and  mid-blade  regions  and  also  at  all  the  rotor 
inlet  nodes.  By  including  all  rotor  nodes,  the  average 
error  with  respect  to  experimental  data  was  27.5%. 

Fig  17  shows  a  plot  of  convergence  criteria, 6,  versus 
the  number  of  iterations  for  a  relaxation  factor  of  0.2. 
The  stability  of  convergence  is  shown  to  initially  decrease 
and  then  after  the  third  iteration  oscillates  about  an 
approximate  value  of  28%.  It  is  important  to  note  that  this 
curve  represents  the  maximum  value  of  6  is  shown  in  equation 
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(II.C.2.5).  In  addition,  the  curve  in  actuality  represents 
the  oscillation  of  nodal  stream  function  values  in  the 
rotor/stator  regions  since  in  fact  this  is  where  the 
non-linearity  is  the  greatest. 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 


Agreement  with  both  experimental  data  and  numerical 
solution  of  Ref . 7  was  best  in  the  stator  region  to  within 
8%.  Predicted  axial  velocity  profiles  in  the  rotor  inlet 
area  were  within  26. 2 J*  of  experimental  results.  The 
instabilities  with  respect  to  static  density  solutions  are 
prevalent.  One  of  the  reasons  for  this  numerical 
disagreement  with  Hirsch  and  Warzee  is  the  isentropic 
assumption  imposed  by  the  present  program.  Recommendations 
for  further  study  on  the  project  include  the  addition  of 
entropy  variations  in  the  rotor  and  stator  blade  regions. 
This  would  necessitate  the  use  of  blade  correlation  data 
[Ref. 13]  for  loss  predictions  and  involve  additional  input 
data  plus  program  additions  to  Subroutine's  SLINE  and  FCAL. 
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APPENDIX   B 
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APPENDIX  C 
SAMPLE   OUTPUT   LISTING 
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APPENDIX  D 


CALCULATION  OF  ROTOR  ELEMENT  FLOW  ANGLES 


The  following  is  a  brief  synopsis  of  the  procedure 
contained  in  Ref.13  for  calculating  the  outlet  relative 
flow  angles  in  a  rotor  element  from  the  given  inlet  relative 
flow  angle  and  blade  solidity.  The  reader  is  referred  to 
Ref.13,  Chapter  VI,  for  specific  details  of  low  speed 
correlation  data. 

As  stated  in  Section  III. A,  uniform  flow  conditions  at 
the  rotor  blade  edges  were  assumed.  This  assumption  coupled 
with  knowledge  of  the  mass  flow  rate  and  rotational  speed, 
enables  one  to  calculate  the  inlet  relative  flow  angle, 3  , 
as  shown  in  Fig  18. 

From  blade  geometry  information,  the  blade  solidity,  5", 

r«     *  (1) 

is  obtained.  At  this  point,  £>  ,  and  T*  are  given  and  one  may 
calculate  B^  ,  the  rotor  outlet  relative  flow  angle  from 
correlation  curves  depicted  in  Ref  13.  The  equation  used  to 
determine  fa,    is  the  following, 


^  =  <z   +  £ 


(2) 


where  K   is  the  angle  between  the  tangent  to  the  Diade   mean 
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camber  line  and  the  axial  direction  (Fig  18)  .  This  is 
obtained  from  the  blade  geometry  data,  <*  is  the  low  speed 
deviation  angle  which  is  obtained  from  the  correlation 
curves  in  Ref.  13.  The  following  equations  show  the 
relationship  between  d  and  the  correlation  data. 


$=     io     +Wa£> 


(3) 


C*(tfL0«Uf. 


ShU^t^  ho  C» 


The  variables  m,K$)  ,k£),  and  S» )  ,  are  all  values 
which  are  obtained  from  the  correlation  curves  and  are  all 
functions  of  the  given  blada  geometry.  The  quantity,  <|>  ,  is 
the  blade  camber  angle  and  again  is  obtained  from  the  blade 
geometry  data.  Once  all  the  variables  are  obtained  from  the 
correlation  data,  equation  (4)  is  solved  for  the  deviation 
angle  for  an  uncambered  blade  section, ^0,  and  then  equation 
(3)  is  solved  for  the  deviation  angle,  ^  .  One  now 
calculates  Sz  from  equation  (2)  for  the  blade  element.  With 
Ax  now  a  known  quantity,  one  now  calculates  the  absolute 
flow  angle,  o(2   ,  from  uniform  flow  assumptions. 

An  example  follows  for  node  numbers  43  and  57  (Fig  6). 
From  Ref.  [12],  Table  II,  the  following  quantities  are 
obtained  assuming  the  angle  of  incidence,  i,  (Fig  13)  is 
zero  and  therefore  the  inlet  relative  flow  angle,  8{  ,  is 
equal  to  k  . 


h 

= 

K 

=   61. 

o 

.38 

r 

- 

i. 

3062 

4> 

= 

6. 

95* 
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£a  =    54.93° 
—  )    max   =    0.035 
tip    radius    =    18.25    in 
hub    radius    =    9.125    in 


Assuming-       uniform        flow      at      the      rotor      inlet      and      a 

rotational    speed   of  4359.5    RPM,    the    following    quantities    are 

determined    from    the    rotor   inlet    velocity   diagram    (Fig    13) . 

_m__  {\ol.  L  lUi/sec)  (lM-4-  ua*/Ft*) 

Vm      =    j>  a     =    — — — 

T  (o.ofc  \bm/FtJ)  TT  (IS.2S2-  %\1*S)  in* 

Vm   =  246.802  Ft/sec 
where  the  area  A,  is  determined  from  the  hub  and   tip   radii 
and  the  density  is  assumed  to  be  0.08  lbm/cu  ft. 

Now  one  is  ready  to  obtain  the  correlation   data.    From 
Ref.[13],   Fig  162,  with  8,  =  61.88*  and  T*  =  1.3062, 

L  J  =  2.50° 
From  Ref.[13],  Fig  162,  with  &x   =  61.38   and  (T=1.3062, 

m  =  0.235 
From  Ref.[13],  Fig  172,  with  t/c) max  =  0.0350, 

ki)  t  =  0.29 
From   Ref.[13],   page   222,   one  uses  the  following  value  of 
(Ki) sh  for  65-series  blades, 

Kd)  sh  =  1.0 

At  this  point  all  the  necessary  data  has  been   obtained   for 
equations  (3)  and  (4)  , 

i0  =  (1.0)  (0.  29)  (2.  50)  =  0.725 
From  equation  (3) , 

S   =  0.725%  0.235(6.95*)  =  2.36° 
Finally,  equation  (2)  gives  the  desired  value  of  p^  , 


1  19 


f 


a  =  54.  93*  +  2.36*=  57.29° 


At  this  point  the  relative  flow  angle  for  node  57  has  been 
obtained,^  =  57.29  .  These  two  values  of  relative  flow 
angles,  A,  =  61.88*  for  node  43  and  A  =  57.29°  for  node  57, 
are  then  read  in  the  program  as  input  data  for  numerical 
computation. 

This  process  is  repeated  at  each  required  blade   element 
section  for  the  proper  outlet  relative  flow  angle. 
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~T    T         AXIAL  DIRECTION,  Z 


Figure  18  -   NOMENCLATURE  FOR  CASCADE  BLADE 
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